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Statistical arguments, seemingly well-justified in higher dimensions, can also be used to derive 
reasonable estimates of Lyapunov exponents x in lower-dimensional Hamiltonian systems. This 
letter explores the assumptions incorporated into these arguments. The predicted x’s are insensitive 
to most details, but do depend sensitively on the nongeneric form of the auto-correlation function 
characterising the time-dependence of an orbit. This dependence on dynamics implies a fundamental 
limitation to the application of thermodynamic arguments to lower-dimensional systems. 
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I. INTRODUCTION 


Computing Lyapunov exponents, which measure the 
average linear stability of chaotic orbits in at oo limit, 
entails solving a matrix harmonic oscillator equation with 
time-dependent frequencies. The differential-geometric 
analogue is a Jacobi equation (equation of geodesic devi¬ 
ation) on a manifold equipped with the Eisenhart metric 
Q. To the extent that aperiodic chaotic orbits can be 
treated as ‘random,’ this exact equation can be approxi¬ 
mated as a stochastic harmonic oscillator equation with 
‘randomly varying’ frequencies @]. The implementation 
of such a picture, along with more detailed assumptions 
that would seem reasonable in a thermodynamic limit, 
has yielded remarkable success in estimating Lyapunov 
exponents in high-dimensional systems 0. In lower di¬ 
mensions, where some of these assumptions are at least 
questionable, one still derives estimates that are qualita¬ 
tively correct; but while the general trends of the dynam¬ 
ics tend to be captured, many details are often missed 

ill. 

The aim of this letter is two-fold, namely (i) to de¬ 
scribe how and to what extent this general approach can 
be improved upon by relaxing some of the underlying as¬ 
sumptions, and (ii) to determine physically why, in lower 
dimensions, much of the dynamics eludes a statistical de¬ 
scription. The analysis will focus on Hamiltonians of the 
form H = + for two representative po¬ 

tentials, namely the Fermi-Pasta-Ulam (FPU) potential. 
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and a higher-dimensional analogue of the dihedral poten- 
tial. 
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II. RELAXING UNDERLYING ASSUMPTIONS 

The key physical assumption underlying the entire 
analysis is that the exact Jacobi equation can be approxi¬ 
mated by a matrix stochastic oscillator equation in which 
the curvature, which determines the frequencies, varies 
‘randomly’ along a chaotic orbit, i.e., that the ‘true’ dy¬ 
namical equation can be replaced by an approximate sta¬ 
tistical equation. Given this assumption, what remains is 
to specify the statistics of the random process. Basically 
this has involved three additional assumptions 0 : 

(1) An assumption of isotropy, which reduces the D- 
dimensional matrix equation to a one-dimensional oscil¬ 
lator equation. 

(2) An assumption that the curvature fluctuations acting 
at any given instant are Gaussian-distributed. 

(3) An assumption regarding the form of the auto¬ 
correlation function and an associated correlation time 
Tc, which characterise changes in the curvature fluctua¬ 
tions along a chaotic orbit. 

Mathematically, the reduction from a U-dimensional 
matrix equation to a one-dimensional stochastic oscilla¬ 
tor follows from the assumption that the configuration 
space can be approximated as quasi-isotropic, i.e., that 
if the Riemann tensor be viewed as the sum of a Ricci 
curvature K and a non-isotropic Weyl projective tensor, 
the latter may be assumed to vanish. Physically, this is 
equivalent to assuming that there are no preferred direc¬ 
tions in space, so that the matrix equation reduces to 
D identical one-dimensional equations. At least for po¬ 
tentials in which different degrees of freedom enter iden- 
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tically, it might seem that this assumption becomes in¬ 
creasingly justified in higher dimensions. Earlier inves¬ 
tigations found that, for some potentials, the analytic 
estimate Xana converges quite rapidly towards the value 
Xnum computed numerically, but that for other poten¬ 
tials it does not 0. This suggests that quasi-isotropy is 
likely not the principal source of error. 

It remains, therefore, to consider the form of the cur¬ 
vature fluctuations both at a fixed instant and as a func¬ 
tion of time. To the extent that the entire phase space is 
chaotic, the form of the curvature fluctuations at a fixed 
instant can be addressed either analytically, given an as¬ 
sumption of ergodicity, or numerically, by evaluating the 
statistical properties of representative orbit ensembles. 
However, it would appear that dynamical issues involv¬ 
ing the auto-correlation function and the correlation time 
must be addressed numerically. 


A. Approximating curvature fluctuation as 
Gaussian-distributed 


The assumption is that curvature fluctuations in dif¬ 
ferent directions can be approximated as nearly inde¬ 
pendent and, at any given instanb Gaussian-distributed. 
For potentials investigated in [^, this approximation 
was completely unjustified in systems with two degrees 
of freedom (dof), and questionable in slightly higher¬ 
dimensional systems. To determine the accuracy of such 
an approximation, we compare the Gaussian constructed 
from numerical estimates of the first two moments of the 
curvature to the exact curvature distribution function 
N[K], which can be computed analytically for potentials 
in which the radius r can be expressed in terms of the 
Ricci curvature K. 

Given the assumption of a microcanonical distribution, 
moments of K satisfy 
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where D is the dimensionality. If one can transform from 
the coordinates (gi, 92 , to) to (A, di, 6 * 2 , 
where K labels surfaces of constant curvature (density), 
eq. 0 becomes 


isoenergy surfaces, the integration extends over the entire 
range of angles Oi . For a H-dimensional potential, a natu¬ 
ral choice for the 9i are the angles V’i of generalized polar 
coordinates (r, ipi ,..., 'tpo-i), with 0 < '0i, ' 02 , ' 0 D -2 < 

TT and 0 < '4>d-i < Stt, in terms of which the distribution 
of curvatures becomes 

2-n TT TT 

N[K] = J d^JD-i j d0z3_2...y d0i^[F;-F]^ 

0 0 0 

X sin^“^ 01 sin^“^-02... sin'0£)_2. (6) 


The results of the direct numerical evaluation of this 
{D — I)-dimensional integral for the dihedral potential 
are shown in Fig. ^ (thick lines). 

If the shape of the curvature distribution does not de¬ 
viate significantly from a Gaussian, i.e., if the higher mo¬ 
ments of the curvature are non-vanishing but small com¬ 
pared to unity, a better estimate can be obtained through 
the asymptotic cumulant expansion 
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where the dependence of 


9q 

d(K,9) 


and V on the trans¬ 


formed coordinates is omitted for brevity. For compact 


Pm{x) is a polynomial of degree m, where m is the highest 
derivative of Z(x) kept in the expansion |(2)). 71 is the 
skewness coefficient and 72 the kurtosis. 

Fig. □ shows three representations Nk{K) of the cur¬ 
vature distribution for the dihedral potential: the Gaus¬ 
sian approximation, the cumulant expansion, and the ex¬ 
act value computed from the eq. along with the 

difference between them as measured by A{Ni,Nj) = 


J {N,{K) - N,{K)f dK 


Given that all the ap- 
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FIG. 1: The distribution of curvatures, N[K], approximated by a Gaussian Z{K) (thin solid line), cumulant expansion 
Ci 2 {K) (dashed line) and exact (thick solid line), as computed by eq. 0, for the dihedral potential with dimensionality 
D — 3, 4, 5, 6 and energy levels E = 1,3,6. The last column shows the absolute difference between the exact value of N[K] 
and the purely Gaussian approximation Z{K) (squares), exact value and the cumulant expansion C\ 2 {K) (diamonds) and the 
Gaussian approximation and the cumulant expansion (triangles). A function of the form {a/'/D) +6 is 6t through data points 
(dashed lines). 


surface, A represents the relative difference between the 
approximations. It is apparent that, as the dimensional¬ 
ity D is increased, the convergence of the exact distribu¬ 
tion towards a Gaussian distribution is of order 
which is the usual thermodynamical behavior. The cu¬ 
mulant expansion provides slightly more accurate esti¬ 
mates to the curvature distribution than the Gaussian 
approximation, but also converges as 1/ '/D. 


B. Auto-correlation function r(t) = {K{0)K{t)) 

The stochastic oscillator equation can be solved analyt¬ 
ically using a method developed by van Kampen 0 if one 
assumes that the correlation time Tc, defined such that 
the auto-correlation function r{t) = {K{0)K{t)) « 0 
for t > Tc, is short compared with the duration of the 
stochastic process. When the correlation time scaled in 
terms of the the average orbital time at that energy is 


short, typically less than t = 2, which is the case for 
the 3-dimensional dihedral potential at low energies and 
the FPU potential for all dimensions and energies, accu¬ 
rate estimates of the Lyapunov exponents are found (Fig. 
0. It is also evident that, for the dihedral potential, the 
correlation time grows with energy, which renders van 
Kampen’s analytic solution increasingly inaccurate. Be¬ 
cause of the long tail of the auto-correlation function, 
even in higher dimensions the agreement between the 
numerical results and the geometric estimates does not 
improve. Empirically, it appears that the threshold at 
which van Kampen’s method yields accurate estimates 
is Tc = i/2 ^ 1. This is consistent with the fact that 
van Kampen’s solution is the hrst term in a power series 
expansion in Tc, with error oc t/. 

In point of fact, estimating Lyapunov exponents as 
outlined in iii assumes that the stochastic process is 
delta-correlated, i.e., r(t) = Ta‘^6{t), where r represents 
a characteristic time scale. However, that assumption 
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FIG. 2: Auto-correlation function F in nnits of dynamical times for the 3-, 6- and 9-dimensional dihedral potential with energies 
E = 0,4, 8,12, (top three rows, respectively), and 4-, 5- and 6-dimensional FPU potential with energies logj^Q E = 1, 1.6, 2.8, 3.4 
(bottom three rows, respectively), along with the analytic (dashed lines) and numerical estimates (solid lines) of the Lyapunov 
exponents (rightmost column). 


was introduced purely on grounds of computational con¬ 
venience, and does not accurately reflect the shape of the 
actual correlation function. Rather, analysis of data such 
as those exhibited in Fig. |21 for both the dihedral and 
FPU potentials, reveals that the auto-correlation func¬ 
tion is well fit by the functional form 


T{t) = Ae ^‘sincnt, (10) 


where 0 ;“^ is comparable to, and strongly correlates with, 
a typical orbital time scale (Fig. |2J). The damping rate p 
appears to correlate with the rate A at which an initially 
localized ensemble of orbits would evolve towards a mi- 
crocanonical equilibrium @ , which implies that it is com¬ 
parable to, but somewhat smaller than, the value of the 
largest Lyapunov exponent - typically A ~ 0.2 — 0.65%. 


C. Characteristic time scale of the process r 

In earlier work, the characteristic time scale r was de¬ 
fined in a seemingly ad hoc fashion as the geometric mean 
of two time scales which were identified using purely ge¬ 
ometric considerations, i.e., + t^^), where 

Ti represents the typical time between successive conju¬ 
gate points on the manifold, i.e., points where the Jacobi 
field of geodesic deviation vanishes, and T 2 corresponds 
to the length scale on which curvature fluctuations be¬ 
come comparable to the average curvature. In point of 
fact, however, Tg and T 2 are closely related to the physi¬ 
cally relevant time scales entering into the problem: quite 
generally, it appears that, to a fair degree of approxima¬ 
tion, Ti oc and, as such, scales as a characteristic 
orbital time scale, and that T 2 oc p~^ and, as such, scales 
as the damping time or, equivalently, the time scale asso¬ 
ciated with chaotic phase mixing as an ensemble evolves 
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FIG. 3: Top three rows represent the 3-, 6- and 9-dimensional dihedral potential, and bottom three the 4-, 5- and 6-dimensional 
FPU potential. First column: geometric characteristic times r (solid lines), ri (short dashes) and r 2 (long dashes). Second 
column: numerically obtained characteristic times ti (solid) and t (dashed lines), as defined in lllll . Third column: characteristic 
time estimate tr scaled to r. Fourth column: ratio of the damping rate p and the largest Lyapunov exponents x- Fifth column: 
for the dihedral and FPU potentials as a function of dimensionality, the Spearman correlation between r and t (crosses) and 
Ti and ti (circles) (first and fourth panel), r and tr (second and fifth panel) and cj and the dynamical time (third and sixth 
panel). 


towards a time-independent equilibrium @ . If we define function of the form m implies a characteristic time 

these two physical time scales and their geometric mean scale 

as 


t2 = -, t ^ = 2(ti ^ +12 ^), ( 11 ) 

P 

we find a strong correlation between r and t, as well as 
between ti and ti (first and fourth panel in the rightmost 
column of Fig. 0 ). The time t 2 is typically an order of 
magnitude longer than ti, which means that the char¬ 
acteristic time t is determined primarily by the shorter 
time scale ti, the orbital time scale. 

Assuming proportionality constants that are compara¬ 
ble for both scalings the geometric mean would thus in¬ 
volve T oc l/(p-|-w). In point of fact, an auto-correlation 





J r(t)tdt 

0 

OO 

/ r(t)dt 
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( 12 ) 


Although this differs from the preceding formula for t, 
we find an excellent correlation between the time tp and 
the time scale t defined geometrically (second and fifth 
row in the rightmost column of Fig. 0 . 
























































6 


III. DYNAMICS 

Statistical estimates of the largest Lyapunov expo¬ 
nent in lower-dimensional Hamiltonian system entail 
two types of approximations: (1) time-independent - 
the space is approximated as isotropic and the curva¬ 
ture fluctuations as Gaussian-distribution, and (ii) time- 
dependent - stochastic oscillations are approximated as a 
delta-correlated process with a characteristic time scale 
T derived using ad hoc geometric arguments. The time- 
independent approximations reduce the dynamics to a 
one-dimensional stochastic oscillator equation; the time- 
dependent approximations are used to solve that equa¬ 
tion via van Kampen’s method. 

At least for the two potentials considered here, the 
time-independent approximations appear well-justified 
and are unlikely sources of significant error. However, 
numerical computations show that the time-dependent 
approximations are not always justified. The simplifica¬ 
tions which reduce the exact integral equation for K(t) 
to a linear first order differential equation for {K{t)) rely 
entirely on the assumption that the auto-correlation time 
tc is short when compared with the duration of the pro¬ 
cess. As is evident from the slow decay of the auto¬ 
correlation function, this assumption typically fails for 
the dihedral potential, so that van Kampen’s method is 
strictly speaking inapplicable. The FPU potential has 
a more rapidly decaying auto-correlation function with 


a much shorter auto-correlation time. In this case, van 
Kampen’s method is applicable and, as illustrated in Fig. 
121 yields more accurate estimates of the largest Lyapunov 
exponents. That this approach works much better for 
the FPU potential than for the dihedral potential may 
reflect the fact that the FPU system is sparsely coupled 
in the sense that each coordinate ‘interacts’ only with the 
neighbouring ones, while in the dihedral potential every 
coordinate is coupled to every other one. This possibility 
would seem related to the fact that many properties of 
dynamical systems appear to depend not simply on he 
number of degrees of freedom, but on their connectance 
0 , a measure of the extent to which these degrees of 
freedom are directly coupled. 

Justifying the time-independent approximations rein¬ 
forces the idea that chaos in lower-dimensional systems 
can be modeled by a one-dimensional stochastic oscil¬ 
lator. However, only in some cases, when the auto¬ 
correlation function damps sufficiently quickly, can van 
Kampen’s method be employed to estimate solutions to 
this stochastic oscillator equation. 
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